Fast Cosmological Parameter Estimation from Microwave Background Temperature 

and Polarization Power Spectra 
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We improve the algorithm of Kosowsky, Milosavljevic, and Jimenez (2002) for computing power 
spectra of the cosmic microwave background. The present algorithm computes not only the tem- 
perature power spectrum but also the E-mode polarization and the temperature-polarization cross 
power spectra, providing the accuracy required for current cosmological parameter estimation. We 
refine the optimum set of cosmological parameters for computing the power spectra as perturbations 
around a fiducial model, leading to an accuracy better than 0.5% for the temperature power spec- 
trum throughout the region of parameter space within WMAP's first-year 3cr confidence region. This 
accuracy is comparable to the difference between the widely-used CMBFAST code (Seljak and Zal- 
darriaga 1996) and Boltzmann codes. Our algorithm (CMBwarp) makes possible a full exploration 
of the likelihood region for eight cosmological parameters in about one hour on a laptop computer. 
We provide the code to compute power spectra as well as the Markov Chain Monte Carlo algorithm 
for cosmological parameters estimation at http://www.physics.upenn.edu/~raulj/CMBwarp 

PACS numbers: 98.70.V, 98.80.C 



I. INTRODUCTION 



The recent Wilkinson Microwave Anisotropy Probe 
(WMAP) first-year data release has shown that, despite 
advances in our ability to compute the power spectrum of 
microwave background anisotropics and increases in com- 
puter power, cosmological parameter estimation is still 
relatively slow and computationally expensive. For ex- 
ample, for the simplest parameter space of flat, Lambda 
Cold Dark Matter (LCDM) models, the WMAP team 
needed about a week of supercomputer time (32 CPUs 
on an SGI Origin) to determine the two parameters joint 
2-er confidence contours in the parameter space 0. Al- 
lowing a non-flat geometry slows down the analysis sig- 
nificantly and more than quadruples the time required 
to draw the 2-er contours. The main bottleneck in the 
parameter estimation process is the theoretical compu- 
tation of the power spectrum anisotropy for a very large 
number of models, scattered throughout the parameter 
space. This is mostly done using publicly available codes 
such as CMBFAST (0) or CAMB (0). To accelerate 
the calculations, (DASH) 1 introduced a shortcut in the 
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calculation of the angular power spectra, based on ana- 
lytical and semi-analytical approximations. 

Rigorous analysis of future microwave background 
and other complementary data sets require recomput- 
ing the likelihood contours for cosmological parameters 
each time new data becomes available. Additionally for 
each new experiment, not one but many likelihood anal- 
yses are required to explore the impacts of potential 
systematic errors. Not every researcher has access to 
months of supercomputer time for these computations; a 
fast and reasonably accurate method to compute confi- 
dence regions for cosmological parameters is needed. Mo- 
tivated by the enormous computational burden that a 
proper 3 to 4-er exploration of the likelihood hypersur- 
face represents and by the related need to understand 
more clearly the physical parameters that control the 
shape of the power spectrum, Kosowsky, Milosavljevic 
and Jimenez (KMJ) Q presented a new set of "physi- 
cal" parameters formed from the usual cosmological pa- 
rameters which have nearly-orthogonal effects on the mi- 
crowave background temperature power spectrum. Thus, 
with this set of parameters, the efficiency of Monte Carlo 
techniques for evaluating the likelihood region increases 
significantly; these parameters have been used in the 
WMAP analysis and then incorporated into Lewis' 
suite of microwave background analysis codes • Subse- 
quent work in a similar vein constructed a set of normal 
parameters that locally transformed the likelihood in pa- 
rameter space into Gaussian form0. 
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In the new KJM "basis" , it is also possible to obtain an 
efficient approximation of the temperature power spec- 
trum in a region around the fiducial model as a simple 
linear extrapolation from a fiducial model power spec- 
trum. This new technique makes possible to compute 
approximate power spectra, several orders of magnitude 
faster than any previous methods, and was extensively 
used by the WMAP collaboration 6] in the preliminary 
analysis of the data. 

Despite a significant increase in computational speed 
and physical insight afforded by the KMJ parameters, 
some remaining technical issues prevented a direct appli- 
cation of the KMJ power spectrum fits to real data. In 
particular: (a) the parameters and the linear fits were 
only computed for the temperature power spectrum and 
not for the polarization; (b) the fiducial model used by 
KMJ had an optical depth to the last scattering sur- 
face which was too small (r = 0.06) compared with 
the WMAP best-fit model (r = 0.17); (c) for a fidu- 
cial model with higher r some covariance remained be- 
tween the slope of the power spectrum n, the amplitude 
of the fluctuations S and the optical depth parameter 
Z = exp(— 2t) which suggested that the KMJ param- 
eters Z, S and n were not completely orthogonal; (d) 
errors in the power spectrum estimation became larger 
than 1% in some relevant regions of the parameter space. 

Here we present a refined version of the KMJ al- 
gorithm, CMBwarp, that solves these problems. We 
extend the method to include polarization spectra, 
make all parameters nearly orthogonal, use a fidu- 
cial model given by the best fit model to WMAPext 
data0, and obtain the power spectra using polyno- 
mial fits for improved accuracy over linear extrapo- 
lations. Throughout the region of parameter s pac e 
within WMAP's first-year data 3-cr confidence region[24|, 
the accuracy is now comparable to the difference be- 
tween CMBFAST and other independent Boltzmann 
codes. We make CMBwarp publicly available from 
http : /www. physics .upenn. edu/^-raulj /CMBwarp; in 
addition, we also release accompanying Markov Chain 
code that can be used to perform parameter estimation. 



II. COSMOLOGICAL PARAMETERS AND 
PHYSICAL QUANTITIES 

We first review the KMJ choice of physical parame- 
ters and describe the slightly modified set of parameters 
presented here which yield an orthogonal basis for tem- 
perature and polarization spectra. Our model space is 
the standard class of inflation-like cosmological models, 
specified by five parameters determining the background 
homogeneous spacetime (matter density ri ma t, radiation 
density f2 ra d, vacuum energy density f^A, baryon den- 
sity fib, and Hubble parameter h), two parameters de- 
termining the spectrum of primordial scalar perturba- 
tions (scalar amplitude S and power law index n), and 
a single parameter r describing the total optical depth 



since reionization. For the present we postpone addi- 
tional complications such as tensors, massive neutrinos 
or a varying vacuum equation of state. Throughout this 
paper, we use a fiducial cosmological model correspond- 
ing to WMAP's best fit LCDM model with parameters: 
h = 0.73, fl m = 0.27, fl b = 0.045, fl v = 0.69, n = 0.99, 
t = 0.166, S — 0.88. As it will be clear below, the 
choice for S is convention-dependent. In our convention 
S = 0.88 which corresponds to as = 0.929 in the fidu- 
cial model. We will construct fast and accurate approx- 
imations to the microwave background power spectra in 
some region of parameter space surrounding this fiducial 
model. 



A. Physical parameters for the Temperature power 
spectrum 

The five parameters describing the background cosmol- 
ogy induce complex dependencies in the microwave back- 
ground power spectra through multiple physical effects; 
therefore their effects on the power spectra are not or- 
thogonal. A characteristic scale in the power spectrum 
is the angular scale of the first acoustic peak. It is ad- 
vantageous to choose this angular scale as a parameter 
which can be varied independently of any other param- 
eters. This angular scale is in turn determined by the 
ratio of the sound horizon at last scattering (which de- 
termines the physical wavelength of the acoustic waves) 
to the angular diameter distance to the surface of last 
scattering (which determines the apparent angular size 
of this yardstick). KMJ used this quantity as the first of 
the following set of orthogonal physical parameters. The 
analytic theory underlying any such choice of phy sical pa- 
rameters has been worked out in detail (see |9l llXlllUll2| . 
and also [HEI). 



The characteristic angular scale is 



A 



(1) 



where a* is the scale factor at recombination, the sound 
horizon is given by 



r s (a) = 



dx 



1/2 



(2) 



and the angular diameter distance in a standard FRW 
spacetime by 



with 



D A (a)=aH- 1 \n~l\- 1 / 2 S k (r), 



{sinr, Q > 1; 
r, Q = 1; 

sinhr, fi < 1; 



(3) 



(4) 
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see, e.g., [l5|. The other physical parameters employed 
by KMJ are: 



B = 




V = 


n A h 2 , 


K = 


Q>* ^mat 


^rad 


M = 





a: 2 tf iad ) 1/2 h 2 ; (5) 



i3 is the physical baryon density determining the baryon 
driving effect on the acoustic oscillations [1G3 (which 
modifies the relative odd and even acoustic peak heights); 
7Z is the matter-radiation density ratio at recombination; 
V determines the late-time Integrated Sachs- Wolfe effect 
arising from a late vacuum-dominated phase, but other- 
wise varying V while holding the other parameters fixed 
represents a nearly exact degeneracy (sometimes called 
the "geometrical degeneracy"); M couples only to other 
small physical effects and is an approximate degeneracy 
direction. This choice of parameters is not unique, but 
this specific set is particularly useful and easy to inter- 
pret. 

Given values for A, B, V, 1Z, and M, they can be in- 
verted to the corresponding cosmological parameters by 
rewriting the definition of A in terms of B, V, 1Z, Ai, 
and h, then searching in h until the desired value for A 
is obtained; Qb and Q\ then follow immediately, while 
f2 ma t and S7 rac j can be obtained with a few iterations to 
determine a precise value for a*. The other cosmolog- 
ical parameters which affect the microwave background 
power spectrum are the optical depth to reionization and 
the amplitude and power law index of the scalar pertur- 
bation power spectrum. These turn out to be crucial pa- 
rameters in order to make an approximate method like 
the one described here to work. For reionization, KMJ 
used the physical parameter 

Z ee er 2 \ (6) 

the factor by which the microwave background 
anisotropies on small scales are damped due to Compton 
scattering by free electrons after the Universe is reion- 
ized. This damping occurs for all scales smaller than the 
horizon size at reionization. 

The primordial power spectrum of density perturba- 
tions is generally taken to be a power law, P(k) oc k n , 
with n = 1 corresponding to the scale-invariant Harrison- 
Zel'dovich spectrum. Making the approximation that k 
and I have a direct correspondence, the effect of n on the 
microwave background power spectrum can be modeled 
as 

C e (n) = Ci(n ) I- , (7) 

which is a good approximation for power law power spec- 
tra. A first-order deviation from a power law, character- 
ized by a — d\nn/dlnk |16| . can be represented in a 



similar way. Note that since the dependence is exponen- 
tial, a linear extrapolation is never a good approxima- 
tion over the entire interesting range 2 < £ < 3000. The 
choice of £q is somewhat arbitrary: changing £q simply 
gives a different overall normalization. In this work we 
choose £o = 550 as the best pivot point to reproduce the 
CMBFAST power spectra. 

These parameters have nearly orthogonal effects on the 
power spectra, so the effect of varying a particular pa- 
rameter can be considered independently from the other 
parameters, within a reasonably large region of parame- 
ter space. This property can be exploited to create sim- 
ple and efficient approximations to the microwave back- 
ground power spectra, based on a few computed power 
spectra. The KMJ paper chose a fiducial cosmological 
model described by a set of parameters Sq and computed 
its power spectrum C^(so) using available codes. The ef- 
fects of parameter A were modeled as a simple rescaling 
of the angular multipole £ and the effects of S as a over- 
all normalization factor, the effects of n as in equation [7| 
above. Then for each of the remaining parameters Sj, a 
linear fit was used to describe the variation of Ce(&o, Si) 
as Si moved away from its value in the fiducial model. 

This set of numerical fits gives approximations to the 
power spectrum which are remarkably good in some re- 
gion of parameter space surrounding the fiducial model. 
However, the linear extrapolations are not good enough 
to compute the power spectrum to the accuracy of the 
CMBFAST code used to compute the fiducial model and 
the dependencies on the individual parameters. Here we 
improve on a simple linear extrapolation by performing 
polynomial fits (to the 4th order) to the computed power 
spectra dependencies. 

Unfortunately, the KMJ choice of physical parameters 
for Z and n does not work if the fiducial model has a 
high r. Even small deviations from the fiducial model 
result in strong covariances between the two parameters 
(see further discussion in Sec. IV). Fortunately, a simple 
technique breaks the degeneracy between these two pa- 
rameters: we use Z as a multiplicative factor to the whole 
power spectrum parameter for £ > 180. For £ < 180 the 
shape of the spectrum does not depend on Z but only on 
n as follows: 

C t (n) = C e (n ) x (^/180) - 8 ("-™°) (8) 

where the subscript denotes the fiducial value. We 
then fix the normalization for £ < 180 by imposing the 
condition that the spectrum must be continuous at £ = 
180. For the amplitude S we use the same definition as 
in KMJ. 

In the rest of the paper, we refer to the parameters 
A, B, V, 1Z, A4, S, and Z as "physical parameters", 
as opposed to the usual "cosmological parameters" fib, 
^m, ^rad, ^vac, h, z r , and a quadrupole-based normaliza- 
tion. The set of approximation techniques for the phys- 
ical parameters, comprising polynomial fits for A, B, V, 



4 




200 400 600 800 1000 1200 1400 




200 400 600 800 1000 1200 1400 



1 .5 



sr 1.0 - 




FIG. 1: Temperature power spectrum comparison between 
CMBFAST and CMBwarp for 10 4 models extracted randomly 
from a MCMC for WMAP first year data. The solid line is 
for flat models while the dashed line is for universes with 
arbitrary geometry. The top panel shows the percent mean 
difference averaged from I = 2 to I = 1500, while the bot- 
tom panel shows the percent RMS difference. Note that the 
agreement is at the 0.5% level or better. 



1Z, and M, a multiplicative Z dependence, and an ap- 
proximate n dependence given by Eqs. Q and J^J pro- 
vide an extremely simple, fast, and accurate method of 
computing the temperature power spectrum of the mi- 
crowave background. Slight additional modifications at 
low i values described below provide excellent fits to the 
E-polarization and the TE cross-power. All of these ap- 
proximations taken as a whole comprise the CMBwarp 
algorithm. 



FIG. 2: Same as fig.^but for the TE power spectrum. Since 
TE crosses zero at several points, the vertical axis is plotted in 
fj,K 2 . The dotted line in the bottom panel shows the cosmic 
variance error for the fiducial model. 



the CMB photons, producing the so called "reionization 
bump" at large angular scales. We model this effect us- 
ing the following fitting formula for I < 16, inspired by 



mg 1 

w 



°' ~ Ce ' l-Z 05 [\ogZ ) (2/100) 

(9) 

where Cf, Eiid is the TE Cy value of the fiducial model 
for each evaluated at 



B. Modeling the polarization power spectra 

For the cross power spectrum between the temper- 
ature and the E-mode polarization, no modification is 
done for £ < 180 and the whole spectrum is modeled 
using the the same prescription for the temperature for 
I > 180 but instead of the multiplicative factor Z we 
use the corresponding polynomial fit for Z. This models 
accurately the TE power spectrum neglecting secondary 
effects. However, reionization of the neutral intergalac- 
tic medium by the first generation of stars re-scatters 



£' = int[(0.105(Z/Z ) 8 + 0.93)4 (10) 

The subscript denotes the fiducial model and int(x) 
denotes the integer part of x; we smoothly interpolate in 
the transition region. This approximation is accurate to 
a few percent. 

For the E-polarization power spectra, the power spec- 
trum for all ^'s is modeled using the same prescription 
as for the temperature power spectrum for £ > 180 but 
instead of the multiplicative factor Z we use the corre- 
sponding polynomial fit for Z. 
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FIG. 3: Same as fig □ but for the EE power spectrum. For I < 100 we show the residual comparison (as in fig. |2J since the 
power spectrum is close to zero for few multipoles. The dotted line shows the cosmic variance error for the fiducial model. Note 
that the agreement with CMBFAST is at the 1-2% level and always below the cosmic variance. 



III. PERFORMANCE OF THE 
APPROXIMATIONS 

The perturbative approach of CMBwarp greatly speeds 
up the calculation of CMB power spectra since it requires 
only a few floating point operations per multipole mo- 
ment calculated. Computation time for a single non-flat 
LCDM model is on the order of 0. 005 seconds for CMB- 
warp; about 100 seconds for CMBFAST on an Athlon 
2400+ workstation. In other words, CMBwarp is about 
20 thousand times faster than CMBFAST. We have com- 
puted our numerical fits up to I of 1500 and tested thor- 
oughly the performance of the fit in the same £ range. 
We have checked that the same type of approximation 
also holds up to £ of 3500. Here we present a direct ap- 
plication to data using WMAP first year data and finer 
angular scale experiments which probe scales up to I of 
1500. 

For CMBwarp to be useful in the analysis of present 
and forthcoming high-precision data sets, it must be not 
only fast but also accurate. Two key questions are: (a) 
can the approximate approach reproduce, with reason- 



able accuracy, CMBFAST power spectra? and (b) can it 
recover unbiased and accurate cosmological parameters 
and confidence regions in a likelihood analysis for a real- 
istic data set? We demonstrate below that the answer is 
yes. 

A. Comparison between CMBFAST and CMBwarp 

The first test we perform is to compare power spectra 
from CMBFAST and CMBwarp for a large set of models. 
For the comparison set, we select models from converged 
Markov Chain Monte Carlos for the WMAP one-year 
data for both flat and non-flat LCDM models. The chains 
are run using CMBFAST (version 4.4 using the high pre- 
cision option) and the Cg are saved; we select models with 
parameters within (roughly) the 3cr marginalized confi- 
dence level region for each parameter. We then calculate 
power spectra for CMBwarp and compare the two sets of 
Cg. For TT, TE and EE power spectra we present both 
the mean difference and the root-mean-square between 
CMBFAST and CMBwarp, as a function of multipole I. 
The first quantity is an estimate of the possible bias intro- 
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FIG. 4: Marginalized 1-D likelihood contours for cosmological parameters recovered from the 1-year WMAP TT dataset with 
CMBFAST(dashed line) and those recovered with CMBwarp (solid line). It is apparent that there is no bias in the recovered 
parameters with CMBwarp. Even for r, which is virtually unconstrained, the recovered CMBwarp likelihood shape is very 
similar to that of CMBFAST. The dotted line is the 1-ct level. 



duced by the perturbative approximation, while the lat- 
ter is an estimate of the scatter. A comparison between 
CMBFAST and Boltzmann codes |l8| for one cosmologi- 
cal model (the standard LCDM model) reveals differences 
at the 0.1% level in the temperature power spectrum and 
the 1% level in the polarization power spectrum. 

Figure shows this comparison for 10 4 flat (solid line) 
and 10 4 open (dashed line) models for the temperature 
power spectrum. The average difference is typically 0.3%: 
of the same order of magnitude of the difference between 
CMBFAST and Boltzmann codes. The RMS plot shows 
not only that CMBwarp is unbiased, but that the RMS 
scatter is only at the 0.5% level. As we demonstrate be- 
low, this level of accuracy is sufficient for accurate recov- 
ery of cosmological parameters from current and future 
CMB temperature experiments. 

The TE power spectrum crosses zero at several val- 
ues of £, so a fractional comparison is not possible. In- 
stead we compute the residual difference in ^K 2 be- 
tween the two models. To be more specific we com- 
pute V(E(^(CMBFAST) - C,(CMBwarp)) 2 )/AU dcls 
(RMS) and £(C £ (CMBFAST) -C e (CMBwarp) )/N modehs 
(bias). Figure [21 shows the results. The discrepancy 



is largest at large I, and is of the order of a 1-2% per- 
cent. These discrepancies are always small compared to 

the cosmic variance ^{Cj T Cf E + {Cf E ) 2 )/(2£ + 1), as 
shown in the figure. 

Finally, fig. 01 shows the comparison between CMB- 
FAST and CMBwarp for the EE power spectrum. For 
10 < I < 100 the EE power spectrum is very close to 
zero, thus for i < 100 we do a residual comparison as 
for TE. The discrepancy with CMBFAST is only few % 
of a [iK 2 and significantly below cosmic variance. For 
i > 100 the EE power spectrum is not close to zero and 
we show a fractional comparison. The agreement with 
CMBFAST is at the 1-2% level, similar to that between 
CMBFAST and Boltzmann codes. 



B. CMBwarp MCMC and parameter estimation 
for WMAP first-year data 

We now turn our attention to the determination of cos- 
mological parameters and their errors. Given the satis- 
factory comparison between our approximate power spec- 
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FIG. 5: Marginalized 1-D likelihood contours for cosmological parameters recovered from the 1-year WMAP TT+TE and 
ACBAR TT datasets with CMBFAST(dashed line) and those recovered with CMBwarp (solid line). Note the good agreement 
between the two likelihoods for all parameters and the absence of a bias in both the best fitting cosmological parameters and 
the likelihood contours. 



tra and those from CMBFAST, we need to demonstrate 
that when we determine cosmological parameters from 
real measurements, CMBwarp yields the same parame- 
ter values and confidence intervals as CMBFAST. 

For this purpose we have calculated likelihood regions 
for two data sets: first the WMAP one-year temperature 
power spectrum data 0, second the WMAP one-year 
temperature and TE power spectra [2C| complemented 
by the ACBAR data points [21j extending to I ~ 1500. 
The above tests are performed for flat and open cosmolo- 
gies. In the likelihood analysis we use the full covariance 
matrix [6| and we explore the likelihood surface with a 
Markov Chain (see |(| for more details). The procedure 
is effectively equivalent to that illustrated in j(| , but with 
CMBwarp replacing CMBFAST. A Markov Chain with 
10 5 models is now computed in 1.5 hours on a single Intel 
3.2 GHz CPU, the computational bottleneck in the cal- 
culation now being the evaluation of the likelihoods and 
not the computation of the power spectra. 

Figure^ shows the marginalized one-dimensional like- 
lihoods obtained from chains computed with CMBFAST 
(dashed line) and CMBwarp (solid line), for WMAP lyr 
TT data only and for a flat model. The recovered best- 



fit parameters display no significant bias (the parameter 
shift is always less that 0.25er), the l-er error is recovered 
to within 15% accuracy or better, and the 2-cr error is 
recovered to 20%. 

Note that despite the fact that t is effectively uncon- 
strained when only temperature data are used, CMBwarp 
still recovers the same shape of the likelihood for r. 

In the second test including the small-scale tempera- 
ture measurements and the WMAP one-year TE data, 
the extra polarization data help to greatly constrain the 
integrated optical depth to reionization. Once again, 
chains of fiat models generated with CMBFAST and 
CMBwarp return virtually the same cosmological param- 
eters as can be seen in Fig. [5j Finally, Fig. [5] compares 
the recovered values of r and h when the model is not 
constrained to be flat. 



IV. DISCUSSION AND CONCLUSIONS 

The set of approximations described here, CMBwarp, 
comprise a highly accurate and extremely fast method of 
approximating the microwave background temperature 
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FIG. 6: Non-flat model 1-d likelihood for h and r using 
WMAP TT +TE. The dashed line is for CMBFAST while 
the solid line is for CMBwarp. 

and polarization power spectra in standard inflation-like 
cosmological models for the range of cosmological param- 
eters which encompass the likelihood regions allowed by 
the WMAP first-year data. For 8 parameters, CMBwarp 
computes each multipole moment with ~ 50 floating- 
point operations; it is unlikely that any significantly more 
efficient computational scheme can be constructed. Note 
that the computational effort is the same for all multi- 
pole moments independent of £; conventional codes re- 
quire large increases in computational cost as larger-^ 
multipoles are computed. This fact will become more 



important in the coming years as experiments such as 
Planck, ACT ^ and the South Pole Telescope map the 
microwave background at high resolution. 

CMBwarp is particularly well-suited to explorations of 
cosmological parameter space using Markov Chain tech- 
niques. As cosmological data continues to improve, joint 
analysis combining many different data sets will continue. 
The technique allows easy incorporation of the true mi- 
crowave background constraints into joint analysis with a 
minimum of computational expense. A useful extension 
of CMBwarp would be a comparably fast and accurate 
computation of the transfer function and matter power 
spectrum. The optimal set of "physical" parameters for 
describing the matter power spectrum will not necessar- 
ily be the same as for the microwave background; the 
crucial property of a set of parameters is that they must 
have orthogonal effects on the matter power spectrum. 
This will be addressed elsewhere. The speed and flexibil- 
ity of CMBwarp make possible analysis with many more 
parameters which are not computationally feasible using 
standard techniques: one example is estimation of the 
primordial power spectrum u sing a number of parame- 
ters to describe it (see, e.g., |23), taking into account 
the effects of the standard cosmological parameters. 

Markov Chain Monte Carlo explorations of cosmolog- 
ical parameter space have now become common. The 
need for very fast, accurate, and portable microwave 
background power spectrum computations in this con- 
text is readily apparent. We hope CMBwarp proves to 
be a useful tool for these analyses. 
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